##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(ArchR)
library(ggsci)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined_peak.combineRNA.qc.Rdata"

    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac_v2"

}


###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
rna_file <- opt$rna_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)
dir.create(paste0(out_path , "/tmp/germ") , recursive = T)
dir.create(paste0(out_path , "/tmp/somatic") , recursive = T)
dir.create(paste0(out_path , "/all") , recursive = T)
dir.create(paste0(out_path , "/germ") , recursive = T)
dir.create(paste0(out_path , "/somatic") , recursive = T)

###########################################################################################

a <- load(comine_data_file)
## testis_combined_peak_combineRNA

b <- load(rna_file)
## scrnat

###########################################################################################
## 细胞颜色
use_colors <- c(pal_npg("nrc")(10) , pal_jco("default")(6))
#names(use_colors) <- unique(scrnat$cell_type)
names(use_colors) <-c( "Myoid cells","Leydig cells" ,          
"Endothelial cells","Zygotene",         
"Round&ElongateS.tids","Patchytene",
"SSC","Sperm" ,                 
"Diplotene","Early stage of spermatids", 
"Leptotene","Sertoli cells",            
"Macrophages","Differenting&Differented SPG",
"Pericytes","NKT cells")

###########################################################################################
## 生殖细胞
cell_level_germ <- c("SSC","Differenting&Differented SPG", "Leptotene",
    "Zygotene","Patchytene","Diplotene",
    "Early stage of spermatids","Round&ElongateS.tids","Sperm"
    )

## 体细胞
cell_level_somatic <- c( "Myoid cells","Leydig cells" , "Endothelial cells", 
"Sertoli cells" , "Macrophages" , 
"Pericytes" , "NKT cells")

###########################################################################################
#### 对于atac 生殖细胞和体细胞重新聚类
## germ
subCells <- getCellNames(testis_combined_peak_combineRNA)[as.character(testis_combined_peak_combineRNA@cellColData[["cell_type"]]) %in% cell_level_germ]
projHeme5 <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = subCells,
  outputDirectory = paste0(out_path , "/tmp/germ"),
  dropCells = TRUE,
  force = TRUE
)

## 重聚类
projHeme5 <- addUMAP( projHeme5 , seed = 1 , force = TRUE)

## atac的聚类图,标记细胞类型
p <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "cell_type", embedding = "UMAP" , pal = use_colors)
p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
image_name <- paste0( out_path , "/tmp/germ/plotEmbedding.atac.pdf")
pdf(image_name , width = 10, height = 8)
print(p)
print(p1)
dev.off()

## 离群的细胞
cell_name <- rownames(projHeme5@cellColData[which(projHeme5@cellColData$Clusters %in% paste0( "C" , seq(1,8) )),])

## 保留的胚系的细胞
use_cell_germ <- subCells[!(subCells %in% cell_name)]

## 输出到对应目录
projHeme5 <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = use_cell_germ,
  outputDirectory = paste0(out_path , "/germ"),
  dropCells = TRUE,
  force = TRUE
)

projHeme5 <- addUMAP( projHeme5 , seed = 1 , force = TRUE)
p <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "cell_type", embedding = "UMAP" , pal = use_colors)
p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
image_name <- paste0( out_path , "/germ/plotEmbedding.atac.pdf")
pdf(image_name , width = 10, height = 8)
print(p)
print(p1)
dev.off()

## 生殖细胞最后保留的细胞
projHeme5_germ <- projHeme5
scrnat_germ <- subset( scrnat , cell %in% gsub( "#" , "_" , use_cell_germ ) )

###########################################################################################
## somatic
subCells <- getCellNames(testis_combined_peak_combineRNA)[as.character(testis_combined_peak_combineRNA@cellColData[["cell_type"]]) %in% cell_level_somatic]
projHeme5 <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = subCells,
  outputDirectory = paste0(out_path , "/tmp/somatic"),
  dropCells = TRUE,
  force = TRUE
)

## 重聚类
projHeme5 <- addUMAP( projHeme5 , seed = 1 , force = TRUE)

## atac的聚类图,标记细胞类型
p <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "cell_type", embedding = "UMAP" , pal = use_colors)
p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
image_name <- paste0( out_path , "/tmp/somatic/plotEmbedding.atac.pdf")
pdf(image_name , width = 10, height = 8)
print(p)
print(p1)
dev.off()

## 离群的细胞
cell_name <- rownames(projHeme5@cellColData[which(projHeme5@cellColData$Clusters %in% paste0( "C" , seq(9,14) )),])

## 保留的胚系的细胞
use_cell_somatic <- subCells[!(subCells %in% cell_name)]

## 输出到对应目录
projHeme5 <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = use_cell_somatic,
  outputDirectory = paste0(out_path , "/somatic"),
  dropCells = TRUE,
  force = TRUE
)

projHeme5 <- addUMAP( projHeme5 , seed = 1 , force = TRUE)
p <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "cell_type", embedding = "UMAP" , pal = use_colors)
p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
image_name <- paste0( out_path , "/somatic/plotEmbedding.atac.pdf")
pdf(image_name , width = 10, height = 8)
print(p)
print(p1)
dev.off()

## 生殖细胞最后保留的细胞
projHeme5_somatic <- projHeme5
scrnat_somatic <- subset( scrnat , cell %in% gsub( "#" , "_" , use_cell_somatic ) )

###########################################################################################
## 所有的
subCells <- c(use_cell_germ , use_cell_somatic)

projHeme5 <- subsetArchRProject(
  ArchRProj = testis_combined_peak_combineRNA,
  cells = subCells,
  outputDirectory = paste0(out_path , "/all"),
  dropCells = TRUE,
  force = TRUE
)

## 重聚类
projHeme5 <- addUMAP( projHeme5 , seed = 1 , force = TRUE)

## atac的聚类图,标记细胞类型
p <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "cell_type", embedding = "UMAP" , pal = use_colors)
p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
image_name <- paste0( out_path , "/all/plotEmbedding.atac.pdf")
pdf(image_name , width = 10, height = 8)
print(p)
print(p1)
dev.off()

## 所有最后用到的细胞
projHeme5_all <- projHeme5
scrnat_all <- subset( scrnat , cell %in% gsub( "#" , "_" , subCells ) )

###########################################################################################
## 输出
## 胚系的
scrnat <- scrnat_germ
testis_combined_peak_combineRNA <- projHeme5_germ
image_name <- paste0( out_path , "/germ/testis_combined_peak.combineRNA.qc.Rdata")
save(  testis_combined_peak_combineRNA , file = image_name )
image_name <- paste0( out_path , "/germ/testis_combined.annotationCellType.qc.Rdata")
save(  scrnat , file = image_name )

## 体细胞的
scrnat <- scrnat_somatic
testis_combined_peak_combineRNA <- projHeme5_somatic
image_name <- paste0( out_path , "/somatic/testis_combined_peak.combineRNA.qc.Rdata")
save(  testis_combined_peak_combineRNA , file = image_name )
image_name <- paste0( out_path , "/somatic/testis_combined.annotationCellType.qc.Rdata")
save(  scrnat , file = image_name )

## 所有的
scrnat <- scrnat_all
testis_combined_peak_combineRNA <- projHeme5_all
image_name <- paste0( out_path , "/all/testis_combined_peak.combineRNA.qc.Rdata")
save(  testis_combined_peak_combineRNA , file = image_name )
image_name <- paste0( out_path , "/all/testis_combined.annotationCellType.qc.Rdata")
save(  scrnat , file = image_name )
